# code to produce weighted regression table and predicted prob plot

# load data

iat_data <- read_rds("./data_clean/rep_file_iat_data.RDS")


##
## descriptive stats
##



##
## regressions
##

iat_data <-
  iat_data %>%
  mutate(city = forcats::fct_relevel(city, "Kharkiv"))


# svy object for weighted regressions
iat_svy_obj <- svydesign(ids = ~ 0,
                         weights = ~new_wt,
                         data = iat_data %>%
                           filter(!is.na(new_wt)))

rus_imp <-
  IAT ~ ethnicity_Russian +
  city + unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3

rus_imp_eth4 <-
  IAT ~ ethnicity +
  city + unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3

rus_exp <-
  explicit ~ ethnicity_Russian +
  city + unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3 +
  IAT

rus_exp_eth4 <-
  explicit ~ ethnicity +
  city + unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3 +
  IAT

# interaction models Kyiv or East City by 4-category ethnicity

rus_imp_eth4_int <-
  IAT ~ ethnicity * city_2_cat +
  unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3

rus_exp_eth4_int <-
  explicit ~ ethnicity * city_2_cat +
  unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3 +
  IAT

# interaction models 2014 parl vote by 4-category ethnicity

rus_imp_eth4_party_int <-
  IAT ~ ethnicity * parl2014_3 +
  city_2_cat +unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3

rus_exp_eth4_party_int <-
  explicit ~ ethnicity * parl2014_3 +
  city_2_cat +unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3 +
  IAT

# interaction models 4 cities by 4-category ethnicity

rus_imp_eth4_int_city <-
  IAT ~ ethnicity * city +
  unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3

rus_exp_eth4_int_city <-
  explicit ~ ethnicity * city +
  unemployed + family_in_Russia + russian_speaking_d + female + age_sd + parl2014_3 +
  IAT

apr_bias_models <- list(rus_imp, rus_imp_eth4,
                        rus_exp, rus_exp_eth4, 
                        rus_imp_eth4_int, rus_exp_eth4_int,
                        rus_imp_eth4_party_int, rus_exp_eth4_party_int,
                        rus_imp_eth4_int_city, rus_exp_eth4_int_city)

# a function that drops variables from models
drop_var <- function(form, vars_to_drop = "- parl2014_3") {
  new_f <- update.formula(old = form,
                          new = str_c(as.character(form)[2],
                                      "~",
                                      as.character(form)[3],
                                      vars_to_drop))
  
  new_f
} 

apr_bias_models[[1]] <- drop_var(apr_bias_models[[1]])
apr_bias_models[[2]] <- drop_var(apr_bias_models[[2]])
apr_bias_models[[5]] <- drop_var(apr_bias_models[[5]])
apr_bias_models[[7]] <- drop_var(apr_bias_models[[7]], "- parl2014_3 - ethnicity:parl2014_3")
apr_bias_models[[9]] <- drop_var(apr_bias_models[[9]])

###


apr_bias_regs <- lapply(apr_bias_models,
                        function(x) svyglm(formula = x, design = iat_svy_obj))

apr_bias_robust_vcov <- lapply(apr_bias_regs, vcovHC, type = "HC")
apr_bias_robust_se <- lapply(apr_bias_robust_vcov,
                             function(x) sqrt(diag(x)))


##
## NEEDED FOR FINAL SEP 04 
##


int_note_new <- c("Robust SEs in parentheses. * p$< 0.05$. ** measured in standard deviation units.")

int_city4_rus_covar_labels_nodem_new <- c(
  "Russian (dummy)",
  "Eth. \\textit{(ref Ukrainian (only))}: Other",
  "\\quad : Ukr. and Rus.",
  "\\quad : Russian (only)",
  "City \\textit{(ref Kharkiv)} : Kyiv",
  "\\quad : Kherson",
  "\\quad : Odesa",
  # "Unemployed",
  #"Family in Russia"
  "Home lang. Russian (dummy)",
  # "Female",
  #"Age**",
  "2014 Vote \\textit{(ref Anti-Russia)} : Pro-Russia",
  "\\quad : Abstain",
  "IAT",
  "Other. eth. $\\times$ Kyiv",
  "Ukr. and Rus. $\\times$ Kyiv",
  "Russian (only). $\\times$ Kyiv",        
  "Other eth. $\\times$ City: Kherson",
  "Ukr. and Rus. $\\times$ Kherson",
  "Russian (only) $\\times$ Kherson",
  "Other eth. $\\times$ Odesa",            
  "Ukr. and Rus. $\\times$ Odesa",
  "Russian (only) $\\times$  Odesa",
  "Constant")



note_all_new <- 
  c("Robust SEs in parentheses. * p$< 0.05$.",
    "Demographics include  \\texttt{family in Russia} (negative in models 3, 4, and 6), \\texttt{age}, and \\texttt{female}.",
    "Observations weighted by raking to voting group.")


tbl_file_name <- "./tables/rus_regressions_all_no_parl_wt.tex"
bias_tex_label <- "tbl:reg_tbl_all_no_parl_wt"

##
## NEEDED FOR FINAL SEP 04 
##

# generates Table E.7, Table F.9

DEMOVAR <- c("family.*|unemployed|female|age.*")
stargazer(
  apr_bias_regs[c(1:4, 9:10)],
  header = FALSE,
  label = bias_tex_label,
  table.placement = "H",
  se = apr_bias_robust_se[c(1:4, 9:10)],
  # type = "text",
  out = tbl_file_name,
  no.space = TRUE,
  # column.labels = rus_column_labels,
  dep.var.labels.include = TRUE,
  dep.var.caption = "Dependent variables: explicit attititudes or IAT $d$-score",
  dep.var.labels = c("IAT", "EXP", 
                     "IAT", "EXP"),
  column.sep.width = "-5pt",
  omit = DEMOVAR,
  omit.labels = rep("Demographics", length(DEMOVAR)),
  omit.yes.no = c("YES", "NO"),
  covariate.labels = int_city4_rus_covar_labels_nodem_new,
  notes = note_all_new,
  notes.align = "l",
  font.size = "scriptsize",
  star.cutoffs = c(.05, NA, NA),
  digits = 2,
  omit.stat = c("f"),
  notes.append = FALSE
)




## preds start here

# SIM CF version

set.seed(20190901)

plot_explicit_sd_units <- TRUE

iat_data_dummy <-
  iat_data %>%
  mutate(other = case_when(ethnicity == "other" ~ 1,
                           TRUE ~ 0),
         Ukr_and_Rus  = case_when(ethnicity == "Ukr and Rus" ~ 1,
                                  TRUE ~ 0),
         Russ_only= case_when(ethnicity == "Russian only" ~ 1,
                              TRUE ~ 0),
         Kherson = case_when(city == "Kherson" ~ 1,
                             TRUE ~ 0),
         Odesa  = case_when(city == "Odesa" ~ 1,
                            TRUE ~ 0),
         Kharkiv= case_when(city == "Kharkiv" ~ 1,
                            TRUE ~ 0),
         proRuss_vote = case_when(parl2014_3 == "Pro-Russia Voters" ~ 1,
                                  TRUE ~ 0),
         abstain = case_when(parl2014_3 == "Abstainers" ~ 1,
                             TRUE ~ 0))

##
## MODEL 9 (implicit bias as outcome variable)
##

m9_formula <- IAT ~ unemployed + family_in_Russia + russian_speaking_d + female +
  age_sd + proRuss_vote + abstain +
  other + Ukr_and_Rus + Russ_only + Kherson + Odesa + Kharkiv +
  other * Kherson +
  other * Odesa +
  other * Kharkiv +
  Ukr_and_Rus * Kherson +
  Ukr_and_Rus * Odesa +
  Ukr_and_Rus * Kharkiv +
  Russ_only * Kherson +
  Russ_only * Odesa +
  Russ_only * Kharkiv

m9_formula <- update.formula(old = m9_formula,
                             new = str_c(as.character(m9_formula)[2],
                                         "~",
                                         as.character(m9_formula)[3],
                                         "- proRuss_vote - abstain"))


# extract data
m9_data <- extractdata(m9_formula, iat_data_dummy)

# fit model, extract PEs and robust VCOV matrix

m9_dum_fit <- lm(m9_formula, data = m9_data, weights = iat_data_dummy$new_wt)
m9_PEs <- m9_dum_fit$coefficients
m9_hcvoc <- vcovHC(m9_dum_fit)

# simulate coefficients
sims <- 10000
m9_simbetas <- mvrnorm(sims, m9_PEs, m9_hcvoc)

##
## MODEL 10 (EXPLICIT bias as outcome variable)
##

explicit_sd <- sd(iat_data$explicit)


m10_formula <- explicit ~ unemployed + family_in_Russia + russian_speaking_d + female +
  age_sd + proRuss_vote + abstain +
  other + Ukr_and_Rus + Russ_only + Kherson + Odesa + Kharkiv +
  other * Kherson +
  other * Odesa +
  other * Kharkiv +
  Ukr_and_Rus * Kherson +
  Ukr_and_Rus * Odesa +
  Ukr_and_Rus * Kharkiv +
  Russ_only * Kherson +
  Russ_only * Odesa +
  Russ_only * Kharkiv +
  IAT

# extract data
m10_data <- extractdata(m10_formula, iat_data_dummy)

# fit model, extract PEs and robust VCOV matrix
m10_dum_fit <- lm(m10_formula, data = m10_data, weights = iat_data_dummy$new_wt) # THIS CHANGED
m10_PEs <- m10_dum_fit$coefficients
m10_hcvoc <- vcovHC(m10_dum_fit)

# simulate coefficients
sims <- 10000
m10_simbetas <- mvrnorm(sims, m10_PEs, m10_hcvoc)

##
##  M9 AND M10 EXPECTED VALUES FOR ALL CITIES, ETHNICITIES, PARL VOTE CATS
##

# create counterfactual scenarios for expected values
all_scens <- 48

other <- c(rep(1, 12), rep(0, 36))
ukr_rus <- c(rep(0, 12), rep(1, 12), rep(0, 24))
rus_only <- c(rep(0, 24), rep(1, 12), rep(0, 12))

kherson <- c(rep(c(rep(1, 3), rep(0, 9)), 4))
odesa <- c(rep(c(rep(0, 3), rep(1, 3), rep(0, 6)), 4))
kharkiv <- c(rep(c(rep(0, 6), rep(1, 3), rep(0, 3)), 4))

abstain <- c(rep(c(1, 0, 0), 16))
pro_russ <- c(rep(c(0, 1, 0), 16))

m9_scens_all <- cfMake(formula = m9_formula, data = m9_data, nscen = all_scens)

# change values of counterfactual scenarios
for(i in 1:all_scens) {
  m9_scens_all <- cfChange(m9_scens_all, "Kherson", x = kherson[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Odesa", x = odesa[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Kharkiv", x = kharkiv[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "other", x = other[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Ukr_and_Rus", x = ukr_rus[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Russ_only", x = rus_only[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "abstain", x = abstain[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "proRuss_vote", x = pro_russ[i], scen = i)
}

# generate predicted values for counterfactual scenarios
m9_preds_all <- bind_cols(linearsimev(m9_scens_all, m9_simbetas, ci = 0.95))

# wrangle data to work w ggplot
colnames(m9_preds_all) <- c("fit", "lwr", "upr")
m9_all_data_ev_plot <-
  bind_cols(m9_scens_all$x,
            m9_preds_all) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               other == 1 ~ "other",
                               Ukr_and_Rus == 1 ~ "Ukr and Rus",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         parl2014_3 = case_when(proRuss_vote == 1 ~ "Pro-Russia Voters",
                                abstain == 1 ~ "Abstainers",
                                TRUE ~ "Anti-Russia Voters"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "implicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)

m10_scens_all <- cfMake(formula = m10_formula, data = m10_data, nscen = all_scens)

# change values of counterfactual scenarios
for(i in 1:all_scens) {
  m10_scens_all <- cfChange(m10_scens_all, "Kherson", x = kherson[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Odesa", x = odesa[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Kharkiv", x = kharkiv[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "other", x = other[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Ukr_and_Rus", x = ukr_rus[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Russ_only", x = rus_only[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "abstain", x = abstain[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "proRuss_vote", x = pro_russ[i], scen = i)
}

# generate predicted values for counterfactual scenarios
m10_preds_all <- bind_cols(linearsimev(m10_scens_all, m10_simbetas, ci = 0.95))

if(plot_explicit_sd_units) {
  m10_preds_all <- m10_preds_all / explicit_sd
}

# wrangle data to work w ggplot
colnames(m10_preds_all) <- c("fit", "lwr", "upr")
m10_all_data_ev_plot <-
  bind_cols(m10_scens_all$x,
            m10_preds_all) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               other == 1 ~ "other",
                               Ukr_and_Rus == 1 ~ "Ukr and Rus",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         parl2014_3 = case_when(proRuss_vote == 1 ~ "Pro-Russia Voters",
                                abstain == 1 ~ "Abstainers",
                                TRUE ~ "Anti-Russia Voters"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "explicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)

m9_m10_ev_all <-
  m9_all_data_ev_plot %>%
  bind_rows(m10_all_data_ev_plot)

m9_m10_ev_all$city <-
  factor(m9_m10_ev_all$city,
         levels = c("Kyiv",
                    "Kherson",
                    "Odesa",
                    "Kharkiv"))

m9_m10_ev_all$parl2014_3 <-
  factor(m9_m10_ev_all$parl2014_3,
         levels = c("Anti-Russia Voters",
                    "Pro-Russia Voters",
                    "Abstainers"))

m9_m10_ev_all$ethnicity <-
  factor(m9_m10_ev_all$ethnicity,
         levels = c("Russian only",
                    "Ukr and Rus",
                    "other",
                    "Ukrainian only"))

ev_full_plot <-
  ggplot(data = m9_m10_ev_all) +
  geom_pointrange(aes(x = ethnicity, y = fit, ymin = lwr, ymax = upr,
                      # linetype = bias,
                      color = bias),
                  position = position_dodge(width = -.5)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed", size = .2) +
  scale_color_manual(values = c("dark gray", "black")) +
  # scale_y_continuous(breaks = seq(-1.5, 1.5, .5),
  #                    limits = c(-1.7, 1.7)) +
  facet_grid(parl2014_3 ~ city) +
  labs(#title = "",
    x = NULL,
    y = "predicted bias (values greater than zero are pro-Ukraine)") +
  theme_tufte() +
  theme(legend.position = "bottom",
        legend.title = element_blank()) +
  coord_flip()

ggsave(filename = ev_fd_rus_plot_name,
       arrangeGrob(exp_ev, imp_ev, exp_fd, imp_fd, ncol = 2),
       width = 7, height = 7)

ev_all_plot_name <- "./plots/ev_plot_all_no_parl_wt.pdf"

# generate Fig F.5
ggsave(ev_full_plot,
       filename = ev_all_plot_name,
       width = 7, height = 7)


